0001    //
0002    //  BigUInt GCD.swift
0003    //  BigInt
0004    //
0005    //  Created by Károly Lőrentey on 2016-01-03.
0006    //  Copyright © 2016 Károly Lőrentey. All rights reserved.
0007    //
0008    
0009    import Foundation
0010    
0011    extension BigUInt {
0012        //MARK: Greatest Common Divisor
0013        
0014        /// Returns the greatest common divisor of `a` and `b`.
0015        ///
0016        /// - Complexity: O(count^2) where count = max(a.count, b.count)
0017        @warn_unused_result
0018        public static func gcd(a: BigUInt, _ b: BigUInt) -> BigUInt {
0019            // This is Stein's algorithm: https://en.wikipedia.org/wiki/Binary_GCD_algorithm        if a.isZero || b.isZero { return BigUInt() }
0021    
0022            let az = a.trailingZeroes
0023            let bz = b.trailingZeroes
0024            let twos = min(az, bz)
0025    
0026            var (x, y) = (a >> az, b >> bz)
0027            if x < y { swap(&x, &y) }
0028    
0029            while !x.isZero {
0030                x >>= x.trailingZeroes
0031                if x < y { swap(&x, &y) }
0032                x -= y
0033            }
0034            return y << twos
0035        }
0036    
0037        /// Returns the [multiplicative inverse of this integer in modulo `modulus` arithmetic][inverse],
0038        /// or `nil` if there is no such number.
0039        /// 
0040        /// [inverse]: https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm#Modular_integers
0041        ///
0042        /// - Returns: If `gcd(self, modulus) == 1`, the value returned is an integer `a < modulus` such that `(a * self) % modulus == 1`. If `self` and `modulus` aren't coprime, the return value is `nil`.
0043        /// - Complexity: O(count^3)
0044        @warn_unused_result
0045        public func inverse(modulus: BigUInt) -> BigUInt? {
0046            var t1 = BigInt()
0047            var t2 = BigInt()
0048            var r1 = modulus
0049            var r2 = self
0050            while !r2.isZero {
0051                let quotient = r1 / r2
0052                (t1, t2) = (t2, t1 - BigInt(quotient) * t2)
0053                (r1, r2) = (r2, r1 - quotient * r2)
0054            }
0055            if r1 >  { return nil }
0056            if t1.negative { return modulus - t1.abs }
0057            return t1.abs
0058        }
0059    }